program PIXGRO

USE mainvariables
USE constants

IMPLICIT NONE

      CALL initconstants  ! initialize

      intrv=3600. !

      droughtflag=0
      irrigationflag=0
      cutcnt=0
      agb=0.
      grain=0.
      harvest=0

      xpsi_1=-50.  ! cmWS
      xpsi_2=-50.  ! cmWS
      xpsi_3=-50.  ! cmWS
      xpsi1=-50.                    ! cmWS
      xpsi2=-50.                    ! cmWS
      xpsi3=-50.                    ! cmWS


      ca=350.  ! what is it?

      finterc=0.1         ! f-Interception
      gapf=0.1            ! gap-fraction



!breed
!breed   emergedate=120
         emergedate=150


!        calculate height dependent air pressure
!        Pressure decreases with height in the first 100 kilometers
!        above the earth's surface according to the formula
!        P(z)=P(sea level)*exp(-z/H). P(z)= pressure at height z,
!        P(sea level)=sea level pressure(~1013 millibars),
!        z=height in meters, H=scale height (to keep the formula simple,
!        we're using 7 kilometers for the scale height)


         atmpress=1013.*exp(-440./7000.)    !440m for center rice

!---------------------------------------------------------------------------------------------------------
!     hrice.for 2006   only the center rice pixel           V 0.2                              02.03.09 gy
!     based on /vol4/korea/UTM/haean/model/h.for from 15.12.08
!---------------------------------------------------------------------------------------------------------
!     18.11.08 sites in korea for UTM:
!             wave rotated   wave not_rotated         fortran      location    bodart   height
!     1  DBF    136,224           136,175             137,176       center       3     459.0000
!     2  DBF     58,317            58,82               59,83         NW          1     778.0000
!     3  grass  122,242           122,157             123,158       center       3     500.0000
!     4  grass  142,122           142,277             143,278        SSE         2     565.0000
!     5  rice   127,205           127,194             128,195       center       3     440.0000
!     6  rice   160,290           160,109             161,110         NE         3     440.0000
!     7  rc     124,256           124,143             125,144       center       3     473.0000
!     8  rc      86,155            86,244              87,245         SW         2     541.0000
!     9  DBF     99,87             99,312             100,313         S          1     754.0000
!
!     todo: remove irrigationflag
!
!---------------------------------------------------------------------------------------------------------
!     intel-compiler8:
!                      ifort -fpe0 -g h1pix.for   (floating underflow results in zero)s
!                      ifort h1pix.for
!---------------------------------------------------------------------------------------------------------


!     read SoilParameter
      call readSoilParas(smatrix,rmatrix)
      call initStartValues(smatrix,rmatrix,1)

!     read CropParameter
      call readCropFiles(cmatrix)

!     read GasFlux-Parameter
      call initseason

!----------------------------------------------------

      open(35,file='hout/p5.asc',status='unknown')    !rice
      open(42,file='hout/t5.asc',status='unknown')    !rice
      write(35,'(a174)')'yy jd hour jd.hh rad tair cantmp pgros darkresp soilresp ctrans trans theta1 theta2 theta3 stemp1 stemp2 &
    &stemp3 lai vint cangl hsum vpd psi1 psi2 psi3 prec lfwt stwt grwt'
!        write(k,'(a174)')'x x x x W/m2 degC degC mmolm-2h-1 mmolm-2h-1
!     + mmolm-2h-1 mm/h mm/h %/100 %/100 %/100 degC degC degC m2/m2 mm/h
!     + mmolm-2s-1 ??? ??? ??? cmWS cmWS cmWS mm/h ??? ??? ???'

      open(2,file='ctrl.dat',status='old')
      open(50,file='met06.asc',status='old') !2006 climate
      read(50,*) header !skip the header line




      do k=1,8760   ! Hour
 20      read(2,*,err=20,end=100)yy,jd,hr   ! read 'ctrl.dat' file
         read(50,*)x,xx,xxx,xdire,xdiff,xtair,xhumid,xwind,xrain  ! read met file

!yy
         droughtflag=0 ! why again??
         cutcnt=0 ! why again?

!        call system('date')
         print*,'***',yy,jd,hr



         xxdire=xdire/1000.     ! W/m2 -> kW/m2
         xxdiff=xdiff/1000.     ! W/m2 -> kW/m2
         humidity=xhumid
         xxtair=xtair
         xxsst=xxtair

         call solar(longitude,rad_lat,jd,int(hr), &
             0.d0,0.d0,cosza_flat,cosza,omega)

         trans=0.
         ctrans=0.
         pgros=0.
         sresp=0.

         bhoe=0.1236*newlai+0.05
         gfac=12.
         gapf=0.1

!        Gasflux crashes if LAI=0.0
         if (newLAI == 0.) newLAI=0.00001

         if (harvest.eq.0) then
!breed      sowdate=90  !inside greenhouse
            sowdate=sstart      !inside greenhouse         !120
         else
            sowdate=harvest+28  !4 weeks after harvest new sowdate
         endif

         call initCropCommon(cmatrix,startww, &
             emergedate,sowdate,laistart)


         if (jd < sowdate.and.harvest == 0) then
            hsum=0.d0
            tsenlai=0.d0
            stage=0
            newLAI=0.00001
            tlai=newLAI
            ww=0.
            lfwt=0.
            stwt=0.
            grwt=0.
         endif
         if (jd == emergedate) then
            newLAI=laistart
            tlai=newLAI
            ww=startww
            lfwt=startww/2.
            stwt=startww/2.
            grwt=0.
            stage=1
         endif

         intercmax=finterc*newLAI
         radsoil=0.

         xxxpsi=(xpsi_2/17.6+xpsi_3)/2.
         yyypsi=xpsi_1
         if (yyypsi.lt.-400.) then !-400.
!ir         3mm irrigation per hour ?
            xrain=xrain+3.
         endif

         call gasflux(rad_lat,jd,xxtair,dble(xwind),humidity,   &
             xxsst,gfac,ca,xxdire,xxdiff,cosza_flat,dble(newLAI), &
            xxxpsi,trans,pgros,atmpress,radsoil,cantmp, &
             darkresp,cangl)

         call nanChecker(jd,hr,pgros,1)

         ctrans=trans

         if (bodart.gt.0) then
            call initSoilCommon(smatrix,rmatrix,bodart)

            call proxel(bodart,intercmax, &
                bhoe,intrv,xxtair,0.1d0,humidity, & !xwind*0.1: wind too high (max=66)
                xxsst,radsoil,dble(newLAI), &
                gapf,atmpress,dble(xrain), &
                trans,sresp,vint,xpsi1,xpsi2,xpsi3, &
                thet1,thet2,thet3,soiltmp1,soiltmp2,soiltmp3)

            call nanChecker(jd,hr,xpsi_3,2)

            xpsi_1=xpsi1
            xpsi_2=xpsi2
            xpsi_3=xpsi3
         else
            sresp=0.
            vint=0.
         endif

         if ((jd == harvest).and.(hr == 0.)) then
            agb=agb+lfwt+stwt+grwt
            grain=grwt
         endif

         call crop(jd,hr,xxtair,pgros,ww,stage,hsum,tsenlai,tlai,pg, &
             rm,dw,newLAI,emergedate,droughtflag,sstart,cutcnt, &
             lfwt,stwt,agb,grwt,harvest,rtwt)

         if (harvest.gt.0.and.newLAI.gt.3.) newLAI=3.
         if (jd.eq.365.and.hr.eq.23.) agb= &
             agb+lfwt+stwt+grwt

         sresp=(((0.32*xxtair+6.)*(-50./((xpsi2+2.* &
             xpsi3)/.9))*exp(190.*(1./(15.+46.)-1./   &
             (xxtair+46.)))))/24.*1000./12. ! g m-2 d-1  -> mmol m-2 h-1

         ewair=6.11*exp(17.26938818*xxtair/(237.3+xxtair))
         ewtp=humidity*ewair/100.
         vpd=ewair-ewtp
         write(35,'(i3,i5,f5.0,f10.5,26f11.3)')yy,jd,hr,   &
             float(jd)+(hr/24.),(xxdire+xxdiff)*1000.,     &
             xxtair,cantmp,pgros*3.6,darkresp*3.6,sresp,   &
             ctrans,trans,thet1,thet2,thet3,soiltmp1,      &
             soiltmp2,soiltmp3,newlai,vint,cangl,   &
             hsum,vpd,xpsi1,xpsi2,xpsi3,xrain,lfwt,stwt,grwt
         write(42,'(i3,i5,f5.0,9f11.4,3i5,f11.4)')yy,jd, &
             hr,newlai,tsenlai,tlai,hsum,lfwt,stwt,grwt,ww, &
             dw,sstart,stage,harvest,rtwt
      enddo

 100  continue
      close(2)
      close(42)
      close(35)
      close(50)

      stop 'Program terminated normally, flush the toilet'


end program PIXGRO
